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An original method to exactly solve the non-Markovian Master Equation describing the interaction 
of a single harmònic oscillator with a quantum environment in the weak coupling limit is reported. 
By using a superoperatorial approach we succeed in deriving the operatorial solution for the density 
matrix of the system. Our method is independent of the physical properties of the environment. We 
show the usefulness of our solution deriving explícit expressions for the dissipative time evolution of 
some observables of physical interest for the system, such as, for example, its mean energy. 
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I. INTRODUCTION 

The problem of the quantum dynamics of a small system interacting with its environment has been extensively 
studied since the origin of quantum mechanics. However, in spite of the noticeable progresses in the theory of open 
I , quantum systems, many conceptual difRculties still remain. Indeed the understanding of the effects of the environment 
' on the physical system of interest, in general, is not an easy task. The conventional way to afford the problem of 
, the description of the time evolution of an open quantum system consists in eliminating the degrees of freedom of 
the environment in order to derive an equation of motion for the reduced density matrix of the small system. This 
\ procedure stems from the fact that we are usually only interested in the dynamics of this subsystem and do not really 
care about the state of the environment and its evolution. The procedure of tracing over the environmental variables 
thus amounts at declaring that everything other than the small system is indifferent to us. In this way one reduces 
the number of variables one needs to take into account to study the problem. 

Only very few physical systems prové to be amenable to an analytic description of their open dynamics. In this 
paper we deal with the simplest open quantum system worth being studied: the damped harmònic oscillator. Such 
a system is of great conceptual importance because it provides a simple, successful starting point in the theoretical 
description of many experimental situations in the quantum òptics, sòlid state physics and quantum field theory 
contexts. For this reason it is also one of the most extensively studied physical systems [Q, g, [| |], g, M. 

Understanding the dissipative behaviour of open quantum systems is a crucial problem for both fundamental 
and applicative reasons. On one hand, indeed, during the last few years, the interpretation of Zurek, associating 
the quantum to classical transition with environment induced decoherence effects, has become increasingly popular. 
• i—i , From this point of view, studying the open system dynamics of exemplary quantum systems, i.e. harmònic oscillator, 
would help in identifying the elusive border between quantum and classical descriptions of the world. On the other 
hand, the huge advances in experimental techniques for controlling the evolution of quantum systems have paved 
the way to the realization of the first quantum lògic gates, the key elements of quantum computers. The biggest 
obstacle to building quantum computers is the decoherence process due to the unavoidable coupling with the external 
environment. Thus, studying a decohering quantum computer is an essential step for the identification of realistic 
quantum error correcting codes Q . 

A Standard method to describe the effects of the environment on a small quantum system is based on the Born- 
Markov Master Equation. Such an approach is vàlid whenever the environmental correlation time is much shorter 
than the typical time seales of the system dynamics. Under these conditions, it is possible to derive and, in some cases, 
to solve analytically the Master Equation ruling the dynamics. This is for example the case of a quantum harmònic 
oscillator coupled through bilincar interactions to a reservoir of harmònic oscillators. This model is a particular case 
of the Quantum Brownian Motion (QBM), where we specify the system to be just a quantum harmònic oscillator ||. 

There exist situations, however, wherein the reservoir correlation time is longer than the system time seales of 
interest and thus the Born-Markov approximation does not hold anymore. This is for instance the case of àtoms 
decaying in photonic band gap materials or atom làsers |ÏC(| . Moreover, very recently, the potential interest of non- 
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Markovian reservoirs for quantum information processing has been demonstrated |ïï{ ] and a non-Markovian description 
of quantum computing, showing the limits of the Markovian approach, has been presented @] . Several procedures 
have been developed to treat non-Markovian processes. Most of them are reviewed in ]ï^| . Very often, in order to 
describe the dynamics of an open quantum system in terms of a master equation for its density matrix, a number 
of approximations involving the relevant system and reservoir parameters are required. One of the most common 
is the weak coupling approximation based on the assumption of a weak system-environment coupling strength. As 
underlined by Paz and Zurek in || , a perturbative approach in the coupling strength can always be shown to lead to 
Master Equations local in time. This fact could appear surprising at a first sight since often non-Markovian quantum 
systems are described by generalized Master Equations involving a non-local memory kernel taking into account the 
past history of the reduced system. The advantage of dealing with a generalized Master Equation which is local in 
time is that the memory effects of the environment are incorporated in its time dependent coefhcients |ï^, [Ï4], fÏ5| . 
This circumstance makes it possible, in some cases, to find an analytic solution of the generalized Master Equation. 

In this paper we present an original method for the derivation of the operatorial density matrix solution of the 
non-Markovian Master Equation for QBM in the weak coupling Hmit. We apply our procedure to the case in which 
the system is a harmònic oscillator interacting with the quantized reservoir through a bilinear coupling. Our method 
does not rely on any assumption other than the weak coupling limit. In particular it is independent of the type of 
environment considered and on the characteristic parameters of the system, provided a perturbative approach to the 
second order in the coupling constant is vàlid. We demonstrate the usefulness and simplicity of our analytic solution 
discussing, as an example, the dissipative dynamics of the mean energy of the harmònic oscillator. 

The paper is organized as follows. In Section II we introduce the mathematical formalism of superoperators and 
we present some key properties for solving the generalized Master Equation for our system. In Section III we describe 
the non-Markovian Master Equation for QBM in the weak coupling limit and in Section IV we derive its operatorial 
solution. In Section V we discuss some approximated forms of the density matrix and we present some applications . 
Finally in section VI we present conclusions. 

II. SUPEROPERATOR FORMALISM 

Let us begin introducing the mathematical formalism and properties which we use in the rest of the paper to 
derive the time dependence of density matrix of the system under scrutiny. Our method is based on the extension of 
the notion of eigensolution of an operator in the 11 superoperator" formalism. 

Given two genèric operators A and p, let us define the superoperators of the S-type A s = (A) s and of the Y,-type 
A s = (^4) s in the following way: 

A s p=[À 7 p], S type (la) 



A^p = {A,p}, Stype (lb) 

where the squared and curi brackcts indicate the commutator and anti-commutator respectively. Eqs. ( p"a| ) and 
( |ïb| ) define a particular class of superoperators hereafter also called commutator (S) and anticommutator (£) su- 
peroperators. It is important to note that a íS-superoperator A s may be linear ly combined with or multiplied by a 
S-superoperator B s giving rise to superoperators belonging neither to the S class nor to the £ class. From Eqs. (p"a|) 



and (lb) linearity follows immediately 

(aÀ + pÈ) = qA s < e ) + [3B S( ^ a,/3eC, (2) 

and the proprieties 

[A S ,B S ] = {À,È} S , (3a) 
[A S ,B S ] = [À,B] S , (3b) 



[A S ,B S ] = [A S ,B S ] = [À,Bf. 



(3c) 
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The main consequence of such relations is that if the commutator of two operators is a c-number then the corre- 
sponding S( E) superoperators do commute.For example, for the super-operators X s ( s ) and P s ^\ corresponding 
respectively to position X and conjugate momentum P operators, the following relation holds 

PCS(S) (P S(E)] = . (4) 

Moreover, from Eq. ( |3c| ) one gets 

[X S ,P S ] = [X S ,P S ] = 2i. (5) 
Another superoperatorial identity extensively used in the rest of the paper is: 

{ÀÈ + ÈÀf = (A S B S + B 5 A S ), (6) 

and, in particular, 

A S A S = (À 2 ) S =A 2S . (7) 

In a certain sense superoperators can be considered as a generalization of the concept of operator used in quantum 
mechanics. Superoperators act on a space whose elements are operators to give other operators. Let us consider a 
generalization of the concept of eigenstates. Looking at the equations 

X S e -,(pX-xP) = _ x e -,(pX-xP) ^ (8a) 

ps ^(px-xp) = _ p e - t (px- x P) j x,p £R, (8b) 

we can identify e l ÍP x ~ xP ) as the eigenoperator of the superoperators X s and V s with eigenvalues — x and — p 
respectively. 

It is not difficult to show that the operators |-^=e l< · p ^ _;ri3 ' ) | define a complete basis (ï|, [ItJ in the operatorial space 

associated to the single particle Hilbert space. Such a basis is orthonormal with respect to the sealar product between 
two genèric operators A and B, defined as follow 

def 



smce 



(A\B) = ' tr{A^B}, (9) 



^-tr{e< p,jt - x ' ps > e -^x-xP) } = S(x - x')S{ P - p'). (10) 



It is then possible to express the density-matrix of a single particle in the following form jï^, [i?], 18 , ^9) 

p = h I x(x,p) e ' l(pÀ ~ xP) dxd p ( n ) 

where the sealar function x( x >P) = tr{e l ( pX ~ xP ^ p} , known in the literature as Quantum Characteristic Function 
(QFC) 0, |ïf, ||, |FJ, satisfies the following properties: 

X(0,0) = 1, x(x,p) = X*(~x,-p), \ X (x,p)\<M (Mel). (12) 

Another superoperatorial eigenvalue equation of interest in this paper is 

N (pX - xP) " = n (pX - xP^j " neN, (13) 
where N = — | (P S X S - X S P S ). In particular we have 

NX n = nX n , NP™ = nP n . (14) 

In the following we will use, for the sake of simplicity, a matrix representation of the previous relations. Let us define 
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With this notation Eqs. (0), ([ÏÏ^) and (|ï^) can be recast in a more compact form 

fs e -x*iÉ = _z e - % *sÈ i (16) 



tr{e lz "' tjÉ e~ l£tJg } = S[z- z% (17) 



2vr 

In these equations J is the following 2x2 matrix 



P=^~ I X(ï)e~ i£tjz d 2 z, (18) 



and thus satisfies the properties J* = J _1 = — J and J 2 = — 1. Note that, using this notation, the superoperator N 
can be cast in the form N = — ~i (Z S )*JZ S . 

III. THE MASTER EQUATION IN THE WEAK COUPLING LÍMIT 

In this section we specify the physical system we wish to study, namely a single harmònic oscillator interacting 
with a quantized environmcnt, and we introduce and discuss the generalized master equation governing its dynamics. 
Let us consider a harmònic oscillator of frequency lüq surrounded by a genèric environmcnt. We express the total 
Hamiltonian H as follows 

H = H + H E + aXÈ, (20) 

where Hq = | [P 2 + X 2 ^j, He and aXÈ are the system, environment and interaction Hamiltonians respectively and 

a is the coupling constant. The interaction Hamiltonian here considered has a simple bilinear form with X position 
operator of the system and E genèric environmental operator. 

For the sake of simplicity we have written the previous expressions in terms of adimensional position and momentum 
operators. Let us denote with p = p the density matrix operator for the oscillator-environment system. 
Let us now assume that 

1. at t = system and environment are uncorrelated, that is p(0) = p(0) (^) /5e(0), with p and ç>e density matrices 
of the system and the environment respectively; 



2. [He, Pe{0)] — (stationarity of the environment ); 

3. trE^E pe(Q)} = (as for example in the case of a thermal reservoir). 

4. a second order perturbative approach in the coupling constant is possible 



Under these conditions one can show |9|, 21 that the non-Markovian generalized Master Equation describing the 

P(t), (21) 



harmònic oscillator dynamics, in the Schròdinger picture, is the following 

dp(t) 



Il·lf -D s (t)+tG s (í) 
in 



dt 

where the superoperator Ds(í) (not S or S type) is defined as 

T> s (t) = í k(t)X s (cosw tX s - sin^oTP 5 ) dr, (22) 
Jo 

Also Gs(t) is neither an S- type nor a S-type superoperator being defined as follows 

G s (í) = / p(r)X s (cos lj tX e - sinw TP S ) dr. (23) 
Jo 
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In Eqs.(|22j) and (J23J) we have introduced the correlation k(t) and the susceptibility /lí(t) functions J22| . Such 
quantities, characterizing the environmental temporal behavior, are defined as follows 

«(r) = ^({È(r), J É(0)}>, (24) 



la 
2ff 



([È(t),È(0)}}, 



(25) 



with (■•■) =tr E {---p E {0)}. 

The form of Eq. ( pï| ) has a clear physical meaning. The superoperator Dg(í), indeed, is strictly related to diffusion 
(decoherence) processes only |J. The superoperator Gs(í), on the other hand, describes dissipation and frequency 
renormalization processes Q . Such a superoperator arises from a quantum mechanical treatment of the environment 
and, indeed, vanishes when the environment is treated as a classical quantity (see also Eq.(p5|)). 

First of all, we concentrate our attention on Gs(í). By using the properties of the superoperators introduced in 
the previous section, one can show that such a superoperator can be recast in the following form 



Gs(t) = \ [r(t)X 2 - 7 (É) (xP + PX)] 5 - t 7 (t) (N + 2) 



where the superoperator N is defined by Eq.(0) and 



r(t) — 2 / /i(r) cos uqt dr, 
Jo 



7(*) 



/Lt(r) sinwoT dr. 



As for Ds(í), it is straightforward to see from Eq. (|22j) that it can be recast in the form 

D s (t) = A(í)(X s ) 2 - n(í)X s P s , 

where 

A(í) = / k(t) cos u)qt dr, 
Jo 

In view of Eq. (||), X s and P s do commute. As a consequence, Dg(t) can be regarded as a quadratic form in the 
commutative operator variables X s and P s . This means that, using the matrix notation introduced in Eq. (|Ï5|), we 
can write 



n(t) = / k(t) sin^oT dr. 
■/o 



(26) 

(27) 
(28) 
(29) 



D s (í) = (Z s )* M(í) Z s , with 



1(í) 



fA(í) 

\-W o . 



(30) 



It is not difficult to check that [Ds(t), Ds(fi)] = 0, whatever t and t\ are. Inserting Eq. ( |26|) into Eq. (21) one 
can write the generalized non-Markovian Master Equation describing the dissipative dynamics of our system in the 
following final form 0, ^ïl 



dp(t) 
dt 



: Hf(í)-D s (í)+ 7 (í) (N + 2) 



P(í), 



with 



ffo(í) = 



À 2 - ^X 2 

CJ 



7(í) 



XP + PX 



(31) 



(32) 



Let us note, first of all, that such a Master Equation is local in time, even if non-Markovian. This feature is typical of 
all the generalized Master Equations derived by using the time-convolutionless projection operator technique ^J] 
or equivalent approaches such as the superoperatorial one sketched in this section. 

Let us have a closer look at the form of our Master Equation. Eq. (^) shows the appearance of two terms, produced 
by the interaction with the environment, modifying the free hamiltonian of the system. The first one, proportional to 
r(t), is a time dependent renormalization of the frequency of the oscillator while the second one, proportional to 7(í), 
describes a coupling between the the X and P operators. The perturbation of the free dynamics due to these terms 
can be easily visualizcd in phase space, as shown in Fig. (Q). Indeed, the effect of the terms proportional to r(t) and 
7 (i) is equivalent to a compression and rotation in phase space of the circle describing the free oscillator dynamics. 
From a dynamical point of view one can show that the term proportional to 7(í) gives rise to both a further frequency 
renormalization and a dynamical dephasing between position and momentum of the oscillator (see Fig. |ï|). Note that 
these features are analogue to the ones present in the dynamics of a classical dissipative oscillator pt] . 
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1: The plain line represents the unperturbed oscillator. The dashed line represents, instead, an oscillator with an hamiltonian of 
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IV. OPERATORIAL SOLUTION 



In quantum mechanics there exists a well established procedure to determine the dynamical evolution of a given 
closed physical system. It prescribes to write down the Schròdinger equation and calculate eigenstates and eigenvalues 
of the Hamiltonian of the system. Once determined these quantities it is possible to describe the dynamics evaluating 
the action of the corresponding temporal evolution operator on the initial state. The dynamics of open systems, 
however, is much more complicated since the description of the state of the system in terms of a state vector is not 
sufRcient anymore. In order to describe completely the physical system one needs indeed to introduce the density 
matrix operator whose time evolution is governed by a master equation definitely more difhcult to handle than the 
Schròdinger equation. 

In this section we present a new method to solve the generalized Master equation, given by Eq. (|3l"|) , for the system 
of interest in the paper. In some sense, the method we describe can be seen as a generalization to open systems of the 
procedure for solving the Schròdinger equation of a closed system, since it is based on the solution of an appropriate 
eigenvalue equation. Indeed we will deal, in the formalism of superoperator, with generalized eigenvalue equations 
involving superoperators and operators, instead of operators and vectors respectively. As we will see in this section, 
the existence of some useful algebraic properties of the superoperators will help us in treating the problem of the open 
system dynamics allowing, in particular, to find the operatorial solution for the density matrix of the system. 

A remarkable virtue of the procedure we are going to describe is its independenec on the expression of the time- 
dependent coefficients appearing in the Master Equation (|3^). This implies that it is applicable to all the Master 
Equations presenting the same structure. 

Our method for solving the Master Equation consists of two steps which can be summarized as follows: 

1) Singling out the temporal evolution superoperator corresponding to Eq.(31); 

2) Understanding how this temporal evolution superoperator acts on p(0). 



A. Temporal Evolution Superoperator 



In the previous section we have introduced the Master Equation related to the dissipative non-Markovian dynamics 
of our system: 



dp(t) 
dt 



Hf(i)-D s (í)+7(*) (N + 2) 



p{t). 



A formal solution of Eq. ( |33| ) can be written as 

m = T(t)p(o) 

where the temporal evolution superoperator T(í) is given by 

" 1 



T(í) = ex Pc - 



^Hf(t 1 )-D s (t 1 )+ 7 (ti)(N + 2) 
in 



dti 



(33) 



(34) 



(35) 
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In (|3^) the pedix "c"stands for the Dyson chronological order. In the following we shall prové that T(í) can be 
factorized as follows: 



T(í) = T 0s (í) T r (t) T D (t), 



where 



To«(í) = exp c 
T r (í) = exp 
T D (í) = exp 



r(í) 



with 



Hf(ti)díj 



(N + 2) 



• / e r(íl) D(Í!) dtí 
Jo 



r(t) = 2 / 7 (íi)díi, 

Jo 

D(í) = Toí(í)D s (t)T .(t). 



(36) 

(37a) 
(37b) 
(37c) 

(38) 
(39) 



To this aim we exploit the Feynman's rule [g3| stating that, whatever the operator or superoperator A(í) and B(t) 
are, 



(A(í) + B(t)) dt 



exp c 



A(t) dt 



exp c 



B(t) dt 



with 



B(t) =[ exp 



A{t) dt 



B(t) exp c 



A(t) dt 



(40) 



(41) 



Applying the Feynman's rule to the time evolution superoperator defined in Eq. ( p5[ ) and taking in consideration 
that, as shown in Appendix A, [N, Hq (í)] = 0, one gets 



T(í)=T Gs (í) T'(í), 



where 



T'(í) = exp c 



j [ 7 (í 1 )(N + 2)-D(í 1 )] dh 
Jo 



(42) 



(43) 



Applying again Feynman's rule in Eq.(f43"l) and using the property (see Appendix A) [N, D(í)l = — 2D(í), leading to 

T r ; 1 (í)D(í)Tr(í)=e r WD(í) I (44) 

we have 



T'(í) =T r (t) T D (t). 



(45) 



Finally, inserting Eq. ( fí5|) into Eq. (^) one obtains the final factorized form of the time evolution superoperator given 
by Eq. @. 

The key advantage of such a factorized form is that it makes easier estimating the action of the time evolution 
superoperator on the initial density matrix of the system. As we shall see in the following, this allows not only to 
solve the Master Equation (|3Ï|), but also to clarify the physical origin of each of its terms. 

It is worth noting that the Dyson chronological order is present in only one of the three superoperators defined 
Eqs. (|36|), (37). Of course, this circumstance leads to a further simplification of the calculations. In fact, while for 
Tr(í) the reason why one can drop the subscript "c"is already clear from Eq. (37), in the case of Td(í) some more 



comment is needed. As we will demonstrate in the following, the Dyson chronological order is not necessary in the 
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expression of Td(í) = exp — e r ^^D(íi) dt\ because the superoperator D(í), defined in Eq.(|39|), is a quadratic 
form in Z s , that is 

D(í) = (Z 5 )* A(í) Z s , (46) 

with A(í) 2x2 matrix of time dependent scalar quantities. Indeed it is straightforward to prové that, when Eq. ( |4^ ) 
holds, [D(i),D(íi)] = 0, for all t,t x . 

In order to prové Eq.(^), let us note that, as demonstrated in Appendix B, the following chain of equalities hold 

T l(t) Z s T Gs (í) = (f l{t)~zf 0s {t)) S = h l{t) Z] S . (47) 



According to these equations, the transformation operated by To s (í) on Z s is a linear transformation correspondent 
to the one operated by To s (t) on Z. One can show (see Appendix B) that 

f l(t) I fos{t) = R(í) É, where R(í) = (<^ t) ^ . (48) 

The functions c(í)and s(t) are the solutions of the Cauchy problems 

2r r(t) 7 2 (í) j(t), „ ,. . fc(0) = l c(0)=0 

y + wl\l- — - - -V }V = with initial conditions ^ ) ( ) ; (49) 

lüq uj^ ^5 ^s(0) = s(0) = w 

The functions s r (í)and c r (í) are defined as 

Cr{t) = m M B{t) _ Sr(t) = m. Ac®. (so) 

Lúç, LÚQ LÜQ LÚQ 

Inserting Eq. ( po|) into Eq. (|3^) and using Eqs. ( fÏ7| ) one gets 

D(t) = (Z s )* M(i) Z 5 , where we have put M(í) = R*(t)M(i)R(i). (51) 
These considerations allow to recast Td(í) in the form 

T D (t) =exp\-{Z s ) t \N{t)Z s ] , (52) 

where 

W(t) = / e^M^díi. (53) 



Jo 

Summarizing, in this subsection we have manipulated the expression of the time evolution superoperator in order 
to put it in a form more convenient for the calculations. In fact we have proved that such a superoperator can be 
factorized in the form given by Eq. (|36|). Moreover we have shown that two of the superoperators appearing in Eq. 
( |36| ) do not need the Dyson chronological order. As we will see in the next section, this circumstance allows to find 
an analytic expression for the density matrix solution. 

B. Evolution of the density matrix 

In this subsection we calculate the action of the factorized superoperator given in Eq. ([36]) on the initial density 
matrix. 



Having in mind Eqs. (p4|) and (36), and writing the initial density matrix as follows (see Eq. (|ÏT|)): 



P(0) = ^ / X o{z')e- l ^ iz d 2 z', ,54 i 



with xa(z') = tr{é z ' tiz p(G)}, we have 



1 



pit) = ^ / T .(t)T r (t)T2,(t) X oC?) e^ z ' iZ d 2 z>. (55) 
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Let us firstly introduce two superoperator eigenvalue equations we will use in the following. Indicating with F(A.) a 
genèric superoperatorial well defined function of the superoperator A, from Eqs. (Q), @ and (^3|) we obtain: 



F(Z £ 



= F(-z), 



F(N)(i*JZ) =F{n){z t }Z 



(56) 
(57) 



where now F(—z) and F(n) are simple scalar functions. 

We begin noting that, using Eqs. (p32|) and (|56|), one can write the term Tb(í) e~ lz ' J ^ appearing in Eq. (^5|) as 



follows 



T D (t) e- lptjf = e~W w(í) P e-"'X 



Moreover, using Eq. (|37|) it is not difficult to prové that 



T r (t) e— ^ = £ [eT( "r JÍ)] " = e^We^H^jl)]. 



(58) 



(59) 



Finally, remembering Eqs. (^) an d ( |4q ) we get 

To s (í) e~ ~ = e ^.w- — e ^, ~ _ e . v -w-v — = e 

Inserting Eqs. (pq) , ( j59f ) and (|60| ) into Eq. (^5|) we can express the density matrix solution as follows 



^ = 2? " 



Xo(z') e 1 (t >exp[e— (-« z'(í) f JZ)] dV 



(60) 



(61) 



After some algebraic manipulation, reported in Appendix D, the previous expression can be recast in the following 
final form 



p(t) = 



2tt 



Xo 



r(t) . 

e — ~R _1 (í) 



2tt 



■2, 



where 



W(í) = e-^ÍR- 1 ®)* W(í) R" 1 ^). 



(62) 



(63) 



Eq. (|62[) constitutes the main result of the paper. It gives the operatorial density matrix solution of the problem of the 
dissipative dynamics of a harmònic oscillator interacting with a genèric environment satisfying properties 2. and 3. 
of section I, as for cxample a thermal reservoir at T temperature. Our approach to the dynamics of the system relies 
on the weak coupling limit but does not invoke the Born-Markov and rotating wave (RW) approximations. In more 
detail, we have solved a non-Markovian generalized Master Equation for the harmònic oscillator coincident with that 
deducible using either the time convolutionless projection operator technique p3[ or the superoperatorial approach 
of pi} ]. At the best of our knowledge this is the first analytical method allowing to solve such a generalized Master 
Equation. In the next section we derive some approximated forms of the solution given by Eq.(|62|) and we show its 
usefulness in calculating the analytic expression of many observables of interest and thus in gaining new insight in 
the dynamics of one of the most extensively studied physical systems. 



APPROXIMATED FORMS OF THE SOLUTION AND APPLICATIONS 



In this section we derive an approximated form of the operatorial density matrix solution, given by Eq. (|62j), vàlid 
when 



Ho(t) « H , 



(64) 



with Ho given by Eq. (|32|). Looking at this equation one sees immediately that Eq. ( p4| ) amounts at neglecting the 
time dependent frequency renormalization and dephasing terms 



r{t)l 2 jt) 



(65) 
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It is worth noting that such an approximation is always justified in the weak coupling regime a <C 1, provided that 
the reservoir frequency cut-off remains finite, as one can appreciate with the help of Eqs. (25), ( p7| ) and (|32]). For 
the same reason it turns out that <C 1. From Eq.(^), it is not difíicult to prové that this last inequality allows 
to write the matrix R(t), defined in Eq. (Eq), as follows 



cos ujot sin uiot 
- s'm ui t coswqÍ 



R(t) 



Inserting this expression in the definition of W(í), given by Eq. (p 



with 



W(í) 



Ca(t) 



,-r(t) 



,r( tl ) 



A(Í X ) , A(*0 



C 2 (í-íi 



(66) 

and exploiting Eqs. ©, (||) and © we get 

n(*i) 



-S 2 (í-íi) 



dtx, 



cos 2woí — sin 2uJot 
- sin 2cjqí — cos 2uot 



Sa(t) 



/ sin 2woí cos 2^0^ 
Icos2üjo^ — sin2woí 



(67) 



(68) 



The form of the matrices defined in the previous equations suggests a further approximation very common in 
Quantum Òptics: the Rotating Wave Approximation (RWA). It basically consists in neglecting rapidly oscillating 
terms as, in our case, those oscillating at frequency 2uiq. This amounts at averaging to zero all the elements of the 
matrices C 2 (í) and S 2 (f) so that Eq. (|67|) reduces to 



W(í) 



,-r(í) 



T(t 



dti 



A r (í) 



2 2 

Substituting Eq. (|69|) into Eq. (|62j) we obtain the following expression for the density matrix solution in the RWA 



(69) 



P(t) 



-RWA 



Zir 



e ^ m Xo 



r(t) , 

e 5-R _1 (í) z 



d'z. 



It is possible to demonstrate that p (t), as given by Eq. (|70|), satisfies the following Master Equation 



r/f 



RWA 



(*) 



A(í) 



fS|2 



+ 7 (í) (N + 2) 



P 



RWA 



(<), 



(70) 



(71) 



which in turn has been derived in pij . 

Once obtained the density matrix we are able, at least in principle, to evaluate the mean value (A) of each and 
every operator A of intercst for the system. One of the advantages of having a solution of the density matrix in terms 
of the characteristic function is the possibility of exploiting the following relations (20[ : 



(X n ) = (-»)» 



(P n ) 



x.p—0 



gn 



x( x ,p) 



(72) 



x.p—0 



By using this relations it is not difíicult to calculate the time evolution of the mean energy of the oscillatory system: 



A r (í) 



with RWA, 



ír{W(í)} without RWA but neglecting renormalization and dephasing terms. 



(73) 



We emphasize that these solutions depend on the initial state of the oscillator only through the term (-ffo)t=o> the 
second term of the sum being independent of the initial state. 

Moreover, looking at Eq. ( |73| ) one sees that the exponential factor accounting for energy dissipation does not 
depend on the initial state of the system but only on the characteristic parameters of the reservoir. 
These features, characterizing Eq. (73), are directly related to the factorized form of the the QCF x*(z). In fact, 
from both Eq. ([32|) and Eq. ([70|) it appears evident that such a function is the product of a gaussian factor and the 

T r(t) "I 

function xo e ~ R _1 (í) z . All the information on the initial state of the system is incorporated in this last function 

from which the first term in the righthand side of Eq. (|7^) comes from. On the contrary, the second term derives 
exclusively from the gaussian factor of the QCF, which depends only on the environment functions and not on the 
initial state of the system. 
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Note from Eqs. (^2|) and ( |70| ) that, whatever the initial state of the system is, 

Xí(z -) ±^ e -(^ W(t) 3 „ (in the RWA) g -^|*1 2 (74) 

This behavior is easily undcrstandable, when the environmcnt is a thermal reservoir, in the light of the thermalization 
process. 

Let us now have a closer look at Eq. (|73|). If we evaluate ír{W(í)} with the help of Eq. @, wc find that 

ír{W(í)} = A r (í). (75) 

This means that the time evolution of the oscillator energy is not affected by the contribution of the rapidly oscillating 
terms neglected in the RWA. This feature comes directly from the particular structure of the frec Hamiltonian operator. 

To better understand this point we consider the gaussian factor appearing in the QCF. The superoperator corre- 
sponding to such a factor has the form 

T G (í) =exp[-(Z s yW(t)Z s ]. (76) 
With some algebraic manipulation one can recast the superoperator appearing in the exponent as follows 

(Z S Y W(t) Z s = [(X s ) 2 + (P s ) 2 ] + ^ [(X s ) 2 - (P s ) 2 ] + 9(i)X s P s , (77) 

where the time dependent coefficients appearing in the previous equation are 

A r (t) = tr{\N(t)}, A(í) = tr{& z \N(t)}, 9(í) = tr{à x W(t)}, (78) 

with <7 Z and à x Pauli spin matrices. 

Exploiting the properties of the trace and putting p'(t) = T(j(í) _1 /3(í), one can easily show that the following chain 
of cqualities holds 

(À) t = tr{p(t)À] = tr{T G (t)p'(t)À} = tr{p\t)T G {t)A] = tr{p> (t)À G (í)} , (79) 
where we have defined 

À G (t)=T G (t)À. (80) 

Having these equations in mind it is not difficult to convince oneself that calculating (A) t using xt(x,p) is equivalent 
at calculating (À G ) t using xo e -- ^ R~ 1 (í) z . In the case À = H n = | (P 2 + X 2 J we have that 



[(X s ) 2 + (P s ) 2 ] H - -2tu , [(X s ) 2 - (P s ) 2 ] H = 0, X s P s ií o = ! (81) 

and thus 

H G (t) = Ho + hujoAr(t). (82) 



This explains why we obtain for 



Note that this equation is not affected by the RWA as well as xo 

the mean value of the oscillator energy the same result with or without the RWA approximation. 

The previous procedure suggests a sufficient condition to single out operators which "do not suffer the RWA 
approximation" , indeed we have that 

[(X s ) 2 -(P s ) 2 ]À = 0, X s P s i = => (À) t = {À)f WA . (83) 

With the expression operators which "do not suffer the RWA approximation" we indicate operators having the property 
that the time evolution of their mean value is not influenced by the RWA. In other words the counter rotating terms 
do not contribute to the dynamics of this class of observables. Examples of operators belonging to such a class are 
X, P, X 2 + P 2 and all linear combinations of such operators, as one can easily verify. On the contrary X 2 and P 2 
do not satisfy condition (p3). For such operators, indeed, exploiting Eq. (|79|), one gets 



(X 2 ) t - (X 2 )f WA = -A(í), {P 2 ) t - (P 2 )f WA = A(t). (84a) 
Similarly, for the mean value of the "correlation" operator (XP + PX) one has 

(XP + PX) t - (XP + PX)f WA =28(t). (84b) 
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VI. CONCLUSIONS 



In this paper we have developed a method to solve the weak coupling generalized Master Equation for QBM. In 
particular we have considered the case in which the system interacting with the environment is a quantum harmònic 
oscillator. The Master Equation we have solved, given by Eq.@, although non-Markovian is local in time. Such a 
Master Equation has been derived in 1 21 by using a superoperatorial technique and in ||, Q by means of the time 
convolutionless method and is the weak coupling approximated form of the exact Master Equation for QBM derived by 
Paz and Zurek in |25|] . The main result of the paper is the derivation of the analytic solution of Eq.(]3Ï|) for the density 
matrix of the reduced system. To this aim we have used an approach based on the àlgebra of superoperators. Our 
method is independent of the specific form of the environment and does not rely on any approximation apart from the 
weak coupling one. We have also studied simpler forms of the density matrix solution obtaincd ncglccting frequency 
renormalization terms and/or performing the RWA. We have demonstrated the existence of a class of superoperators 
whose mean value is not affected by the presence of the counter rotating terms at any time t and we have given a 
sufhcient condition to verify if a given operator belongs to such a class. This circumstance simplifies substantially the 
calculations since, for operators belonging to such a class, one can use the approximated density matrix solution, given 
by Eq.([ro|), in order to calculate their mean value. The analytic solution we derive and discuss in the paper is given 
in terms of the QCF (see Eqs.(|62|)) by means of which one can calculate the expectation value of many observable of 
physical interest in a very direct way, as suggested by Eqs.(|7^). For examplc, thanks to the simplicity of the analytic 
solution we have derived, we succeed in calculating the dissipative time evolution of the mean energy of the system. 
Finally it is worth noting that from the QCF it is easy to derive the Wigner function characterizing the state of the 
dissipative system. 

Concluding, we believe that, the non-Markovian analytic approach we have derived in this paper for QBM can be 
generalized to other fundamental dissipative systems, such as for example the Jaynes-Cummings model with losses. 
As for the QBM, we think that the analytic solution of the density matrix may be used for studying new aspects of 
such a bàsic model, both for fundamental and for applicative research, under conditions in which, up to now, only 
numerical approaches were possible. 
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APPENDIX A: PROPERTIES OF DISSIPATION SUPEROPERATOR 

In this Appendix we demonstrate some useful properties of the superoperator N. Such properties, allowing to 
factorize the time-evolution superoperator given in Eq (pa), are the following 



[N,fi£(f)]=0, [N,D(Í)] = -2D(i). 



(Al) 



The first one is the most surprising one: it can be considered as a "dynamical invariance" of the dissipation process. 
As we shall see in the following, the "dynamical invariance" originatcs from the particular structure of the Hamiltonian 

tcrm H (t), defined by Eq. (||). 

Independently on the form of the time dependent functions present in such a term, it can be written, in general, as 



H (t) = Z l S(t)Z with S*(í) = S(í) 
where S(í) is a 2 x 2 real symmetric matrix. Remembering that 



one gets 



[N,H S (Í)] 



N=-- Í (Z E ) Í JX^, 



(Z S ) 4 JZ S , (l*S(í)Íy 



(A2) 
(A3) 



■-% (Z S )*J 



Z l S{t)Z,Z 



Z l S{t)Z,Z 



s\ 4 



jz- 



(A4a) 
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It is possible to show that [E6j 



Z É S(í)Z,Z =tJS(t)Z. 



(A4b) 



Substituting this equation into Eq. (A4a) and using the properties of the matrices S(í) and J (cfr Eq. (|Ï9|)) onc 
obtains 



[N,Hf(í)] = -i ((Z s )*JJS(í)Z s + (Z s ) t S t (í)J t JZ s ) 



-- ( -{z^ys(t)z b + (z L· )'s(í)z ò ) = o. 



(A4c) 



The generality of the matrix S(í) ensures the validity of this relation for a genèric quadratic hamiltonian and thus, 
in particular, for Ho(t). 

As far as the second commutation relation is conserned, remembering that D(í) = (Z s )* M(í) Z s one has 



[N,D(í)] 



N, (Z s )* M(í) Z' 



N,Z 



M(í) Z ò + {Z b f M(í) 



N,Z f 



Exploiting the superoperatorial commutation rules given in Eq. (0) yields 



N,Z 4 



P S X S -X S P S ) ,Z & 



= z k 



Finally, substituting into Eq. ( |A5a| ) one gets 

[N,D(t)] = -(Z S Y M(t) Z s - {Z s ) f M(í) Z s = -2D(í) 



(A5a) 



(A5b) 



(A5c) 



APPENDIX B: TIME EVOLUTION 

In this Appcndix we shall demonstrate that the general time evolution of a "S"- or "S"-type superoperator is 
equivalent to the "S"- or "S"-type superoperator of the time evolution of the corresponding operator. In formulas 
this amounts at demonstrating Eq. (EtI). First of all let us define T c (t) as solution of the equation 



-T c (t) = L(t)T c (t) 

where in general [L(t), L(t\)] ^ 0. L(t), as well as T c (í), can be either an operator or a superoperator. 
The previous equation can be solved in iterative way and its solution is 



n=0 Jo Jo 



L(íi) • • • L(t n ) dti--- dt n = exp c 



L(ti) dU 



It can be shown that T c 1 (t) satisfies the equation j2^] 

jT-^t) = -T-\t)L{t) 
which can be again solved in iterative way giving the following form of the solution 



n=0 



L(t n ) ■ ■ ■ L(íi) dti ■ ■ ■ dt n = exp a 



L(íi) dti 



(Bl) 



(B2) 



(B3) 



(B4) 



Let us define 



A c (t)=T c (t)AT- 1 (t)=T c (t)A 



(B5) 
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From the previous equation we have 



^-T c (t)A= [L(t),A c (t)] =L s (t)A c (t) =L s (t)T c (t)A 



-T c (t) = L s (t)T c (t) 



(B6) 



This means that T c (t) satisfies an equation similar to the one given by Eq. (Bl). The form of T c {t) can be obtained 
from Eq (|B2|) replacing T c (t) -> %{t) and L(t) -> L s (t). 

Now, let us suppose that A is "S"- or "S"-type superoperator and L is an "S"-type superoperator. From the 
properties given in (|^) we have that 

A s c ^(t) = T e (t) A 5 < s ) T-^í) = T c (t)A s ^ 

= E / •■■ [ n ~ 1 [L s (t 1 )---[L s (t n ),A s W}---}dt 1 ---dt n 



n=0 



n=0 



E 
E 

exp c 



o Jo 



í»-i 



'[L(ii) •••[£(*„), i] ••·] s ^díi·.·dí B 



S(S) 



L s (í 1 )···L s (í„)Adí 1 ···dí„ 



L s (íi) díi 



S(S) 



(Tc(í) i) 



S(E) 



S(S) 



(ic(í)) 



S(S) 



(B7) 



where we have put T c (í) = exp c / Q L s (ti) dt\ 
Now, let us consider 



Af< s >(í) - T- X (í) A 5 ( s > T c (í) = T c 7\t)A s ^. 



(B8) 



The last definition comes directly from the definition of T c (t) and can be verified applying T c (t) on Af(t) and, 
viceversa, appl ying the previous definition of T^ft) on Af (í). 



Exploiting Eq.(B4) and following the same lines of the derivation of Eq. (B7), one can show that 

S(S) , , , s S(S) 



Af( E )(í) = (f- 1 (í)if c (í)) =(i (í))' 



Let us now derive Eq. (|48|). Using Eq. ( A4b ), (Bl), (B3) and remembering the definition of Z(t): 

h)=T3l{t)Èf 0s {t). 



hase 



-Z{t) = --folit) [H (t),Z] T 0s (t) = SS(t)Z(t) S(í) = -S(í). 



Deriving the previous equation once more we get 



^Z(t)= -JS(í) + (JS(í)) 2 Z(t)= -JS(í) - det[S(t)] Z\t). 



In our case 



S(í) 



1 /wo-r(É) 7 (í) 



2 V 7(í) w 

so that the differential equation for the operator X(t) has the following form 

dí^ V w / 



(B9) 



(B10) 



(Bll) 



(B12) 



(B13) 



(B14) 
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A solution of the previous equation can be written as 

X(t) = Xc{t) + Ps(t), 



(B15) 



with c(í) and s(t) solutions of Eq. (49) with the same initial conditions. Note that, from (Bll), it follows that 



-X(t)=u P(t)+>y(t)X(t). 



With th help of Eq. (B15) we have for P(t) 



P(t) = LUxQ-T®*® 



LÚQ 



= Pc r (t) - Xs r (t), 



where we have defined, as in Eq. (50) 



-(*) = 



c(t) - j(t) 3 (t) 

LJq 



S r (t) 



c(í)- 7 (í)s(í) 



LÚQ 



Grouping Eqs. ( B15 ), ( B17 ) and using the matrix representation we obtain 



X(t) = Xc(t)+Ps(t) 
P{t) = -Xs r (t) + Pc r {t) 



Z(t) = R(t), R(í) = 



c(t) s(t) 

-S r (t) C r (t) 



(B16) 



(B17) 



(B18) 



(B19) 



APPENDIX C: FINAL FORM OF THE DENSITY MATRIX 

In this ap pendix we sketch the main steps of the derivation of the final form of the density matrix solution, given 
by Eq. (|62|), from Eq. (|6Ï|). Let us remind the genèric expression for the density matrix in terms of the QCF (See 
Eq. (|Ü|)) 



(Cl) 



where Xt(z) = tr\e lz JZ /5(í)}. From Eq.(|6l"|) we have 



P{t) 



2tt 



-W'w(í) z' n (J\ P r(t) 



Xo(z') é (t >exp[e— (_j z >(ty SZ)] d 2 z' 



Comparing these last two equations one gets 



Xt (z) = tr{e lZti2 p{t)} = ±-j e-^'WCt) e m tr { e * ^z exp [ e ^ (_ t P(t)*j|)]} d 2 z> 

From Eq.@ and using some properties of the delta function we obtain 



-^trie^'^expie^i-izityiÉ)]} = 6[z- e^z'it)] = e~ T{t] 5[z' - e^R" 1 ^) z\. 



Introducing such this equation into Eq. (|C3|) yields 



(C2) 



(C3) 



(C4) 



X t(z) = I w(í) ? xq{z')5[z' - e-^R-^t) z\ d 2 z< 

= expi-e^izf (R- 1 (í)) í W(í)R- 1 (í) z\ X v[e~^ R" 1 ^) z\ 



(C5) 
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Finally, substituting this equation into Eq.(Cl) gives 




/ 



e 



(25* W(í) z 



r r(t) _ 1 

Xo e 2 R x (í)2 e 




(C6) 



where 



W(í) = e-^ÍR-H*))* W(t) R" 1 ^). 



(C7) 



which is the form of /5(í) given in Eq.(|G2|). 
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